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Abstract. With the help of computer algebra we study the diagonal matrix elements (Or p ), where 
O = {l,/3,ian/3} are the standard Dirac matrix operators and the angular brackets denote the 
quantum-mechanical average for the relativistic Coulomb problem. Using Zeilberger's extension of 
Gosper's algorithm and a variant to it, three-term recurrence relations for each of these expectation 
values are derived together with some transformation formulas for the corresponding generalized 
hypergeometric series. In addition, the virial recurrence relations for these integrals are also found 
and proved algorithmically. 



Science is what we understand well enough to explain to a computer. Art is everything else we do. 

cr 

Donald E. Knuth [22] 

1. Introduction 

This work has been initiated by the following email regarding Doron Zeilberger's Z60 conference: 

O 

(N 



^ . http://www.matri.rutgers.edu/events/Z60/ 



Email to Peter Paule from Sergei Suslov [27 Feb 2010] 
Subject: Uranium 91+ ion 

"... I understand that you are coming to Doron's conference in May and write to you with an unusual 
suggestion. . . 

I am attaching two of my recent papers inspired by recent success in checking Quantum Electrodynamics 
in strong fields [see Refs. [36] and [38] in this paper]. 

It is a very complicated problem theoretically, and fantastically, enormously complicated (at the level of 
science fiction!) experimentally, which has been solved - after 20 years of hard work by theorists from 
Russia (Shabaev + 20 coauthors/students) and experimentalists from Germany. 

Experimentally they took a uranium 92 atom, got rid of all but one electrons, and measured the energy 
shifts due to the quantization of the electromagnetic radiation field! 



Date: June 10, 2012. 

1991 Mathematics Subject Classification. Primary 81Q05. Secondary 33C20. 
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Mathematically, among other things, the precise structure of the energy levels of the U 91+ ion requires 
the evaluation of certain relativistic Coulomb integrals, done, in a final form, in my attached papers . . . 

Here is the problem: 

These integrals have numerous recurrence relations found by physicists on the basis of virial theorems. 
They are also sums of 3 (linearly dependent) 3F2 series. 

Now you can imagine what a mess it is if one tries to derive those relations at the level of hypergeometric 
series (3 times 3 = 9 functions usually!). 

It looks as a perfect job for the G-Z algorithm in a realistic (important) classical problem of relativistic 
quantum mechanics. It looks as a good birthday present to Doron, if one could have done that. I feel 
we can do that together. 

Looking forward to your answer on my crazy suggestion, BW, Sergei" 

The first named author's computer algebra response reported at the Z60 conference is presented 
in this joint paper. 

2. Relativistic Coulomb Integrals 

Recent experimental and theoretical progress has renewed interest in quantum electrodynamics 
of atomic hydrogenlike systems (see, for example, [3], [10], [11], [13], [14], [17], [32], [34], [35] and 
the references therein). In the last decade, the two-time Green's function method of deriving formal 
expressions for the energy shift of a bound-state level of high-Z few-electron systems was developed 
[32] and numerical calculations of QED effects in heavy ions were performed with an excellent 
agreement to current experimental data [10], [11], [34]. These advances motivate a detailed study 
of the expectation values of the Dirac matrix operators multiplied by the powers of the radius 
between the bound-state relativistic Coulomb wave functions. Special cases appear in calculations 
of the magnetic dipole hyperfine splitting, the electric quadrupole hyperfme splitting, the anomalous 
Zeeman effect, and the relativistic recoil corrections in hydrogenlike ions (see, for example, [1], [31], 
[33], [36] and the references therein). These expectation values can be used in calculations with 
hydrogenlike wave functions when a high precision is required. For applications of the off-diagonal 
matrix elements, see [23], [24], [25], [29], [30], and [33]. 

Two different forms of the radial wave functions F and G are available (see, for example, [19] and 
[39]). Given a set of parameters a,ai,a 2 , (3, /3 l , /3 2 , and 7, depending on physical constants s,k,/i, 
and v, consider 

\G{r) ) y 7 r(n + 2z/) 1 1 J \ [5 1 [3 2 ) \ L^{2a(5r) ) 1 ' 

where, using the notation from [20], L^(x) stands for the corresponding Laguerre polynomial of 
order n. Throughout this paper, 

K = ±(j + 1/2) , v = Vk 2 -H 2 , 

f x = aZ = Ze 2 /hc, a = Vl-e 2 , (2.2) 
e = E/mc 2 , (3 = mc/h, 
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and 7 = — v){sk — v), with the total angular momentum j = 1/2, 3/2, 5/2, etc. (see [4], [5], [8], 
[26], [36], and [39] regarding the relativistic Coulomb problem). The following identities 

e/i — a (u + n) , eji + au = a{n + 2v) , eji — au = an, (2.3) 

2 2 2 2 / , o \ 2 22 

e k —v— an [n + 2u) — /i — a k 
are useful in the calculation of the matrix elements. 

The relativistic Coulomb integrals of the radial functions, 

POO 

A p = / r p+2 (F 2 (r) + G 2 (r)) dr, (2.4) 
Jo 

POO 

B p = / r p+2 (F 2 (r) - G 2 (r)) dr, (2.5) 

/•oo 

C p = / r p+2 F (r) G (r) rfr, (2.6) 
Jo 

have been evaluated in Refs. [36] and [38] for all admissible integer powers p, in terms of linear 
combinations of special generalized hypergeometric 3_F 2 series related to the Chebyshev polynomials 
of a discrete variable [18], [19]. 

Note. We concentrate on the radial integrals since, for problems involving spherical symmetry, one 
can reduce all expectation values to radial integrals by use of the properties of angular momentum. 

Throughout the paper we use the following abbreviated form of the standard notation of the 
generalized hypergeometric series 3 F 2 ; see, e.g., [20]: 

/ ai, a 2 , a 3 \ f ai, a 2 , a 3 \ ™ MMMk , 



h, h ) ' \ 6i, b 2 ' J ^ (6i)fc(&2)fcA;! 



where (a)* := a(a + 1) ... (a + k — 1) denotes the Pochhammer symbol. 

Analogs of the traditional hypergeometric representations for the integrals are as follows [36]: 

T(2v + l) ( 1 — n, —p, p+1 \ 

2^ (2a(5f -i { + ' . A p = 2pean 3 F 2 ' " P (2.8) 

^ v ^ r(2z/ + P + i) p y y 2i/ + i, 2 y 

1 — n, — p, p + 1 \ / — n, — p, p+1 

+ (p + a K ) 3 F 2 \ 2j/+ij t j +&.-«.) 2j/+ij i 

(1 — n, — p, p + 1 \ / — n, — p, p+1 

2u + l, 1 J y 2u + l, 1 

4n(2aP) p -I^±1L gp (2.10) 



r(2z/ + p + i) p 

2z/ + 1, 1 / v ' ' " " \ 2z/ + 1, 1 



1 — n, — p, p + 1 \ / — n, — p, p+1 

= a (/i + a/c) 3 F 2 - a (// - a/c) 3-F2 
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The averages of r p for the relativistic hydrogen atom, namely the integrals A p , were evaluated in 
the late 1930s by Davis [6] as a sum of certain three ^F 2 functions. 1 But it has been realized only 
recently that these series are, in fact, linearly dependent and related to the Chebyshev polynomials 
of a discrete variable [36]. The most compact forms in terms of only two linearly independent 
generalized hypergeometric series are given in Ref. [38]. 

In addition, the integrals themselves are linearly dependent: 

(2k + e(p + 1)) A p - (2ek + p + 1) B p = 4/iC p ; (2.11) 

see, for example, [1], [29], [30], and [36]. Thus, eliminating, say C p , one can deal with A p and B p 
only. 

The integrals (2.4)-(2.6) satisfy numerous recurrence relations in p, which provide an effective 
way of their evaluation for small p. A set of useful recurrence relations between the relativistic 
matrix elements was derived by Shabaev [30] (see also [1], [7], [29], [33], [36], and [40]) on the basis 
of a hypervirial theorem: 

2kA p -( P +1)B p = 4{iC p + 4f3eC p+1 , (2.12) 
2kB p -{p+\)A p = 4(3C P+1 , (2.13) 
fiB p -(p + l)C p = f3(A p+1 -eB p+l ). (2.14) 

From these relations one can derive (see [1], [30], and [33]) the linear relation (2.11) and the following 
computationally convenient recurrence formulas (2.15)— (2. 18), stated in our notation as 

A - (n I ]) ^ 2 e + 2K(p + 2) + e(p+l)(2Ke + p + 2) 
Ap+1 ~ ~ {P + 1) 4(1 -e 2 ) (p + 2)ftz Ap (2 - 15) 

4^ 2 (p + 2) + (p + 1) (2ks + p+l) {2ks + p + 2) 
+ 4(l- £ 2 )(p + 2)p> p ' 

B , - (y I n ^ 2 + 2^(2p + 3)+5 2 ( P +l)( P + 2) 

Bp+1 ~ {P+1 > 4(1-5^ + 2)/^ Ap (2 - 16) 

4/x 2 £ (p + 2) + (p + 1) (2ke + p + 1) (2k + e [p + 2)) 
+ A(l-e 2 )(p + 2)(3fj J p 

_ A^e (p + 1) + p (2ke + p) (2k + e (p + 1)) 

" ^ /i(4z/ 2 -p 2 )p Ap {2 - U) 



and 



q 4/x 2 (p + 1) + p (2/te + p) (2/te + p + 1) 

^ 71 2 2\ Pi 

li [Av z — p z ) p 



r /3 4z/ 2 + 2/ tg (2p + l) + g 2 p(p+l) 

fi \Av — p z ) 
4u 2 e + 2k (p + 1) + ep (2ke + p + 1) 
^ /x(4z/ 2 -p 2 ) p ' 



respectively. 



1 He finishes his article by saying: "In conclusion I wish to thank Professors H. Bateman, P. S. Epstein, W. V. Hous- 
ton, and J. R. Openheimer for their helpful suggestions." 
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Note, (i) These recurrences are complemented by the symmetries of the integrals A p , B p , and C p 
under the reflections p — > —p — 1 and p — > —p — 3 found in [36]; see also [2]. (ii) These relations 
were also derived in [37] by a different method using relativistic versions of the Kramers-Pasternack 
three-term recurrence relations. 

3. Computer Algebra and Software 

The general algorithmic background of the computer algebra applications in this paper is Zeil- 
berger's path-breaking holonomic systems paper [41]. The examples given in the following sections 
restrict to applications: (i) of Zeilberger's extension [43] of Gosper's algorithm [9], also called 
Zeilberger's "fast algorithm" [42, 22], and (ii) of a variant of it which has been described in the un- 
published manuscript [21]. Both of these algorithms have been implemented in the Fast Zeilberger 
package zb.m which is written in Mathematics and whose functionality is illustrated below. A very 
general framework of Zeilberger's creative telescoping (i), and also of its variant (ii), is provided by 
Schneider's extension of Karr's summation in difference fields [12]; see, for instance, [27, 28] and 
the references therein. 

The Fast Zeilberger Package can be obtained freely from the site 

http : //www . rise . jku . at/research/combinat/sof tware/ 

after sending a password request to the first named author. Put the package zb .m in some directory, 
e.g., /home/mydirectory, open a Mathematica session, and read in the package by 

In[l] := SetDirectory ["/liome/ppaule/RISC_Comb_Sof tware_Sep05 . dir/f astZeil"] ; 
<<zb .m 

Fast Zeilberger Package by Peter Paule and Markus Schorn (enhanced by Axel Riese) 
- © RISC Linz - V 3.53 (02/22/05) 

A Mathematica notebook containing a full account of the Mathematica sessions described below, 
together with some additional material, is available at: 

http : //hahn . la . asu . edu/~ suslov/ curres/index . htm 



4. Unmixed Three-Term Recurrence Relations 

The following relations purely in the A p and B p , respectively, have been established in [38]: 

A 1 = » P (P) A (A1) 

p+1 a 2 (3(^ 2 (p+l)+p(2€K+p)(2en + p+l))(p + 2) p K ' J 

(4z/ 2 - p 2 ) (4/i 2 (p + 2) + (p + 1) {2sk +p+l) (2sk + p + 2))p 



(2a/3) 2 (4/i 2 (p + l)+p (2ek + p) (2en + p + 1)) (p + 2) 

o £»Q(P) o (A 9) 

p+i ~ a 2 P (4z/ 2 + 2sk (2p + 1) + e 2 p (p + 1)) (p + 2) v K ' 



(4z/ 2 - p 2 ) (4z/ 2 + 2en (2p + 3)+e 2 (p+ 1) (p + 2)) (p + 1) 
(2a/3) 2 (4z/ 2 + 2en (2p + 1) + e 2 p (p + 1)) (p + 2) 



B 



p— i) 
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where 

P(p) = 2ep(p + 2)(2en + p)(2eK + p+l) (4.3) 
+e [4 (eV -v 2 )-p (4eV +p (p + 1))] 
+ (2p + 1) [4e 2 K + 2(p + 2) (2e/2 2 - «)] , 

Q(p) = (2p + 3) [4z/ 2 + 2«;(2p+l)+p(p+l)] (4.4) 
-a 2 (2p + 1) (p + 1) (p + 2) . 

In comparison with other papers (e.g., [1], [2], [29], [30], [36], [37], and the references therein), this 
approach provides an alternative way of the recursive evaluation of the special values A p and B p , 
when one deals separately with one of these integrals only. The corresponding initial data A — 1 
and B_i = a 2 (3/ ' ji can be found in [36]. 

Note. The derivation in [38] resembles the reduction (uncoupling) of the first order system of 
differential equations for relativistic radial Coulomb wave functions F and G to the second order 
differential equations; see, for example, [19] and [39]. 

With Zeilberger's definite extension [42, 43] of Gosper's algorithm [9] for indefinite hypergeometric 
summation, the derivation of such recurrences is fully automatic if the input is given as a terminating 
hypergeometric series (and provided that the input is of computationally feasible size). We illustrate 
this by a mechanical derivation of the following simple three-term recurrence relation for the integral 
C p , not found in [38]: 

C p+l - M2P + 1) a 2 (3{p 2_ AK 2 Hp+1) C p ( 4 - 5 ) 

(p2 _ [(p + 1)2 _ 4k 2] 

+ P 9 C v -1. 

(2a/3) 2 (p 2 - 4/t 2 ) (p + 1) P 

As input for C p we take the hypergeometric sum representation from (2.10). We start our Mathe- 
matica session by reading in the RISC "Fast Zeilberger" package: 

In[l] := «zb.m 

Fast Zeilberger Package by Peter Paule and Markus Schorn (enhanced by Axel Riese) 
- © RISC Linz - V 3.53 (02/22/05) 

In[2] := (a_) k _ := Pochhammer [a,k] ; 

F1 T> 1 = (~P) k (P +1 )k F9 r k -, = (~n) k (~p) k (p + l) k 

L J-_ (2 v+l\ (l) k k! ' t2l ^-~ (2,+ l) k (l) k k! ' 
In [3]:= FullSimplify[{Fl[k]/Fl[k] , F2 [k] /Fl [k] }] 



Out [3]= {1, -— — } 
k — n 



Gamma|2 v + II \ 1 r , / . . . , n 



In [4] := f [kj := 4 /x (2 a ^ —\ *Fl[k] * a(// + a/c)-a(/i-aK 



Gamma[2 v + p+l]/ '" \ " n — k 
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In[5]:= SuslovRec=Zb [f [k] , k, 0, Infinity, p, 2] // Simplify 



Out [5] = {4 a (3 (-(3 + 2 p) (-(2 + 3 p + p 2 ) fi 2 (n + v) - 2 a n k \x (n + 2 u)+ 



a 2 k 2 (n + i/) (2 + 4 n 2 + 3 p + p 2 + 8 n v)) SUM[l+p] + 

a (2 + p) (-(1 + p) 2 / u 2 + a 2 k 2 (4 n 2 + (l + p) 2 +8 n v)) SUM[p + 2] == 

(1+p) (l + 2p + p 2 -4i/ 2 ) ( (2 + p) 2 fi 2 + a 2 k 2 (4 n 2 + (2+p) 2 + 8 n v)) SUM[p]} 



Here C p = SUM[p]. Utilizing two of the identities (2.2)-(2.3) brings Out [5] into the form (4.5). In 
order to prove the correctness of Out [5] , just type 

In [6] := Prove [] 

and the program generates automatically a pretty print version of a proof in a separate window or 
file, respectively. 

The computerized derivations and proofs of (4.1)-(4.2) are analogous; one finds the details in the 
corresponding Mathematica notebooks on the article's website. 

5. Related Transformations of Generalized Hypergeometric Series 

Several relations between two pairs of the generalized hypergeometric series under consideration 
are given in [36] and [38]: 




(5.1) 



(2u + n)(2u + p + 1) (2z/ + p + 2) (2w + p + 1) 
4i/ (2z/ + 1) (u + n) (p+1) 




and 




(5.2) 



n (4z/ + 2n - p - 1) (2u + p + 1) (2i/ + p + 2) 
4i/ (2i/ + 1) (y + n) (p + 1) 
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In addition, 



Pip 




2v + n 

_ (p-2i/)(2i/ + p + l) 
2 (2i> + 1) (1/ + n) 3 2 



(5.3) 




and 



p(p+l) 

3-^2 





(5.4) 

These relations are "responsible" for the transformation between two different hypergeometric forms 
of the relativistic Coulomb integral [36, 38]. The second named author was able to give only the 
proof of the last relation from the advanced theory of generalized hypergeometric functions. 

With the zb.m package, one can not only prove but also find such relations, in the literature 
called also contiguous relations, automatically. We illustrate this by a computer derivation of (5.1). 

In[l] := «zb.m 

Fast Zeilberger Package by Peter Paule and Markus Schorn (enhanced by Axel Riese) 
- © RISC Linz - V 3.53 (02/22/05) 

In[2]:= (a_) k := Pochhammer [a,k] ; F0 [k_] := (l ~ n) * ( ~ P)k ^ ± , l)k ; 

(2 v+ l) k (l) k k! 

F1[k] . = (l~n) k (~p-l) k (p + 2) k - p2 _ = (-n) k (-p- l) k (p + 2) k - 
(2 ^ + 2) k (l) k k! " (2 v\ (l) k k! 

Our goal is to compute rational function coefficients Co, ci, C2, free of the summation variable k, such 

that 

oo 

(c F0[k] + Cl Fl[k] + c 2 F2[k}) = 0. 

k=0 

With a parameterized version of Gosper's algorithm, similar to Zeilberger's extension of Gosper's 
algorithm, we compute such q together with a hypergeometric expression g such that for non- 
negative integer N: 

N 



k=0 

This is accomplished using the option Parameterized; finally we send N to infinity: 
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In [3] := FullSimplify[{l, Fl [k] /FO [k] , F2 [k] /FO [k] }] 

(l+k + p) (1 + 2.) nO+k + p) (k + 2,) 

1 (-1 + k-p) (l+k + 2 v)' 2 (k-n) (-1+k-p) v 5 

In [4]:= Gosper [FO [k] , {k, 0, N} , 

- - f< (l+k + p) (1+2 i/) n (1 +k + p) (k + 2 v) ln 

Parameterized -> {l, -- — — -, — ^ . y r^j] 

1 (-1+k-p) (l+k + 2 v)' 2 (k-p) (-1 + k-p) z/ 

If C N' is a natural number, then: 

N 

Out [4]= {J^4 (1+p) v (n + u) (1 + 2 v) F [k]-(l + 2 n + p) (n + 2 v) (l+p + 2 v) 

k=0 

(2 + p + 2 v) Fi[k]+ 2 n i/ (1+2 v) (1 + 2 n + p + 4 i/) F 2 [k] == 

- ((1+N + p) (1 + 2 v) (2n + 4n 2 + nN + 2n 2 N + 3np + 2n 2 p + nNp + 
np 2 +4 v + 8 n v + % n 2 i/ + 2 N i/ + 6 p v + 8 np z/ + 2 Np z/+ 
2 p 2 1/ + 8 v 2 + 8 n v 2 + 8 p i/ 2 ) 

Pochhammer [1 - n, N] Pochhammer [-p, N] Pochhammer[l + p, N])/ 
((l+N + 2 v) N! Pochhammerfl, N] Pochhammer[l + 2 i/, N]) } 

For — > oo this gives the desired relation because the right hand side is when N > p. 

The computerized proofs of (5.2)-(5.4) are similar and the corresponding Mathematica notebooks 
are available on the article's website. 

6. Virial Recurrence Relations 

A general procedure of verification of the linear relations between the relativistic integrals can be 
formulated as follows. Start from the hypergeometric series representations for the integrals involved 
into the identity/relation in question, and find all linear dependencies between the corresponding 
hypergeometric series using the package zb.m. Substitute the integrals into the desired identity, 
eliminate the linear dependent sums/vectors from this equation, and then simplify the coefficients 
in front of the rest of the series to zero with the help of the standard identities among the quantum 
numbers of the relativistic Coulomb problem. 

One can easily see that the linear relation (2.11) is equivalent to (5.4), and that (2.12) follows 
from (2.11) and (2.13). To illustrate our strategy, we derive (2.13) directly from the hypergeomet- 
ric representations for the relativistic Coulomb integrals (2.8)-(2.10). To this end, we input the 
hypergeometric summands involved in the relations (2.8)-(2.10) and (2.13): 

In[l] := «zb.m 
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(l-n) k (-p) k (p + l) k 



In[2]:= (a_) k _ := Pochhammer [a,k] ; F0 [k_] : = 



(2 u+ l) k (2) k k! ' 



FIDO := (l 7 n)k ( - p)k / P + , l)k ; F2[k_] := ^ ™* <* + ^ ; 

(2 i/+l) k (l) k k! (2 i/ + l) k (l) k k! 

F3 r k1 •= (l~n) k (-p-l) k (p + 2) k - p4 _ = (-n) k (-p - l) k (p + 2) k 

(2 i/+l) k (l) k k! ' " (2 !/ + l) k (l) k k! 



In [3] := FullSimplify[{l, Fl [k] /F0 [k] , F2[k]/F0[k], F3[k]/F0[k], F4 [k] /F0 [k] }] 

„ ut[3: . {1 , 1 + k , _<*±iLs, 1+k+ ?A£±i), <L±jE) ;» Wp> } 

k — n 1-k + p (k-n) ( — 1 + k-p) 

r r n r ! r (l +k) II 2 k ( 1 + k) 

In [4]:= Gosper [F0 [k] , {k, 0, N} , Parameterized -> {1, 1 + k, — , 1+kH -, 

k — n 1— k + p 

(1+k) n (1+k + p) 
(k-n) (-l+k-p) 1 

If C N' is a natural number, then: 

N 

Out [4]= {^-n Fi[k] + (l+n + p) F 2 [k]-n F 3 [k] + (-1 + n - p) F 4 [k] ==0, 



k=0 
N 



^2 n p F [k] + (l + 2 n + p + 2 u) F 2 [k] + (-1 - p - 2 v) F 4 [k] == 



k=0 



2 n (l + N + p) Pochhammerfl — n, N] Pochhammerf— p, N] Pochhammerfl + p, N] 
N! Pochhammer[2, N] Pochhammer[l + 2 u, N] 

N 

^2 -2 n (n + 2 v) Fi[k] + (l + 2 n + 2 n 2 + 2 p + 2 n p + p 2 + 2 v + 

k=0 

4 n v + 2 p v) F 2 [k] - (1 + p) (1 + p + 2 v) F 4 [k] == 
2 n (1+N) (l+N + p) Pochhammer[l — n, N] Pochhammer[— p, N] Pochhammer[l + p, N] 
N! Pochhammer[2, N] Pochhammerfl + 2 u, N] 



} 



Notice that for N — > oo all the right hand sides vanish because they are when N > p. 

Summarizing, the package has found the following three linear relations: 

n (X + U) - (1 + n + p) Y + (1 - n + p) V = 0, (6.1) 

2np Z + (1 + 2n + p + 2z/) F - (1 + p + 2u) V = 0, (6.2) 

2n(n + 2v)X + (1 + p)(l + p + 2is)V (6.3) 
= [(n + l) 2 + 2p + {n + p) 2 + 2(2n + p + F 
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for the following five linear dependent vectors 

1 - n, -p, p + 1 \ / —n, -p, p+1 

" - 1 2i/ +1, 1 ; ^ 2i/ + 1, 1 

" • ' 2i/ + 1, 2 J' 3 2 I 2i/ + 1, 1 



— n, — p — 1, p + 2 
2i/ + l, 1 

We choose to present everything in terms of Y and V: 

In [5] := Linl = n*(X + U) - (1 + n + p)*Y + (1 - n + p)*V ; 

Lin2 = (2 n p)*Z + (l + 2n + p + 2 z/)*Y - (1 + p + 2 z/)*V ; 
Lin3 = -2 n (n + 2 z/)*X + 

(l + 2n + 2n 2 + 2p + 2np + p 2 + 2z/ + 4nz/ + 2pz/)*Y 

(1 + p) (1 + p + 2 z/)*V ; 
Solve [Linl==0 && Lin2==0 && Lin3==0, {X, U, Z}] ; 
FullSimplifyE /,] 



Out [51= ffX^ -(1+P) V (l+P + 2 z/)+Y (2 n 2 + 2 n (l+p + 2 i/) + (l + p) (l+p + 2 z/)) 
11 2 n (n + 2 v) 

V (2 n 2 -2 n (l+p-2 z/) + (l + p) (l+p-2 z/))-(l+p) Y (l + p - 2 v) 



U 



2 n (n + 2 i/) 
V (l+p + 2 u)-Y (1 + 2 n + p + 2 i/) 

2 n p " 



Introducing the Coulomb integrals, 

In [6]:= Ap[X, Y, Z] := (2 // (2 a /3)~ (p) (Gamma[2 i/ + 1])/ (Gamma [2 z/ + p + 1]))~(-1)* 
((/i + a k)*X + (/x - a k)*Y + 2 p £ a n*Z) ; 

Bp[X, Y, Z] := (2 fi (2 a /3)~ (p) (Gamma[2 v + 1])/ (Gamma [2 v + p + 1]))~(-1)* 
(e (/i + a k)*X + e (// - a k)*Y + 2 p a n*Z) ; 

CpluslDJ, V] := (4 /x (2 a /3)~(p + 1) (Gamma[2 v + 1])/ 

(Gamma[2 v + p + 2]))~(-l)*(a (// + a k)*U + a (/x - a k)*V) ; 



we express the desired relation in terms of X, . . . , V: 

In [7]:= (2 K)*Bp[X, Y, Z] - (p + l)*Ap[X, Y, Z] - (4 /3)*Cplusl [U, V]; 

7, /. Gamma [2 + p + 2 z/] -> (l+p + 2 z/)*Gamma[l + p + 2 z/] ; 
FullSimplif y [%] 
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Out [7]= - 2~ 1 -P (a 

H Gammafl + 2 v) 

(/i ((X + Y) (1 + p - 2 £ re) + U (1 + p + 2 v) - V (1 + p + 2 v)) + 
a(2npZ(£ + pe-2re)+re((X-Y)(l + p- 2£re) + 
U(l + p + 2z/)+V (l + p + 2 i/)))) Gamma[l + p + 2 i/] 

Next we rewrite the relevant part into a linear combination of X, ... ,V: 

In [8] := ZERO = (/x ((X + Y) (1 + p - 2 e re) + U (1 + p + 2 v) - V (1 + p + 2 u) ) + 
a(2npZ(e + p£-2re) + 

re ((X - Y) (1 + p - 2 e re) + U (1 + p + 2 v) + V (1 + p + 2 v)))) ; 
Collect [ZERO, {X, Y, Z, U, V}] ; 
FullSimplif y [*/.] 

Out [8]= 2anpZ(£ + pe-2re)+Y(l+p-2£re)(-are + / [x) + 

X (1 + p - 2 e re) (a re + /x) + V (a re - /x) (1 + p + 2 z/) + U (a re + /i) (1 + p + 2 i/) 
Eliminating X, U and Z, 

-(1 + p) V (l+p + 2 u)+Y (2 n 2 + 2 n (l + p + 2 i/) + (l+p) (l+p + 2 u)) 



In [9] := '/. /. {X 
U 



2 n (n + 2 i/) 

V (2 n 2 -2 n (l+p-2 z/) + (l + p) (l+p-2 i/))-(l+p) Y (l + p-2 v) 



2 n (n + 2 i/) 
V (l+p + 2 u)-Y (1 + 2 n + p + 2 v) 

O ^ ^ ' 

2 n p 

FullSimplif y[°/„] ; 
Collect ['/., {Y, V}] ; 
FullSimplif y [%] 

(1+p) V (l + p + 2 u) (-u (n-e re + z/) + a (n 2 e-n re + £ re 2 + 2 n £ - re v)) 

Out [9]= — ; v v ^ '-!-+ 

n (n + 2 v) 

Y ( (l+p-2 s re) (-a re + „) - (1 + P) ( » " + (l + P " 2 "> (l+P + 2 "> 



2 n (n + 2 i/) 
a (e + p e-2 re)(l + 2 n + p + 2 v)+ 

(l + p-2 e re) (a re + /x)(2 n 2 + 2 n (l + p + 2 i/) + (l+ p ) (l+p + 2 v)) 

2 n (n + 2 v) 
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Finally, we simplify the coefficients of V and Y: 

In [10] := ZeroV = - / u(n-ere + z/)+a (n 2 £-nre + £re 2 + 2nei/-rei/) ; 
ZeroY = 2n(l + p- 2ere) (-a k + fi) (n + 2 v) - 
(1 + p) (a re + jx) (1 + p - 2 v) (1 + p + 2 v) - 
2 a n (e + p e -2 re) (n + 2 ^) (1 + 2 n + p + 2 v) + 

(1 + p - 2 e re) (a re + /x) (2 n 2 + 2 n (1 + p + 2 v) + (1 + p) (1 + p + 2 u)) ; 
{ZeroV, ZeroY} ; 
FullSimplify[%] ; 
% /. n ->(£//- a i/) /a ; 
FullSimplif y [%] ; 
7, /. e~2 - > 1 - a~2 ; 
FullSimplif y [%] ; 
7, /. re~2 - > v~2 + [i"2 

Out [10]= {0, 0} 

which is the name of the game. 

Computerized proofs of (2.14) and some of its extensions work the same; they are available on 
the article's website. 

In a similar fashion, seeking for a more general linear combination of the corresponding integrals, 
with the zb . m package one can derive the following two-parameter relation: 

[D{p + 1) - C (2k + e(p + 1))] A p - [2D re - C {2en + p + 1)] B p + 4fiC C p + A(5D C p+l = 0, (6.4) 

where C and D are two arbitrary constants. The virial relations (2.11)— (2.13) are its special cases. 

We would like to point out the following relation: 

(p + 1) (2reC7 p - M P ) = Pip + 2) (sA p+1 - B p+1 ) , (6.5) 

as another simple example. 

Note. This relation is a linear combination of (2.12)-(2.14); see [1]. 

7. Conclusion 

The relativistic Coulomb integrals (2.4)-(2.6) were recently evaluated in a hypergeometric form 
[36]. The corresponding system of the first order difference equations (2.15)-(2.16) has been solved 
in [38] in terms of linear combinations of the dual Hahn polynomials thus providing an independent 
proof. Here, with the help of the Fast Zeilberger Package zb.m we give a direct derivation of these 
results. 

One of the goals of this article is to demonstrate the power of symbolic computation for the 
study of relativistic Coulomb integrals. Namely, computer algebra methods related to Zeilberger's 
holonomic systems approach allow not only to verify some already known complicated relations, 
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but also to derive new ones without making enormously time-consuming calculations by hands or 
with ad hoc usage of computer algebra procedures. 

In a sequel to this article we are planning to investigate the computer-assisted derivation of 
recurrences, e.g. the "birthday recurrences" from Section 4, by taking as a starting point the 
original definition of the Coulomb integrals (2.4)-(2.6). To this end, we will use Koutschan's 
package HolonomicFunctions [16]. This package, written in Mathematica, implements further ideas 
related to Zeilberger's holonomic systems paradigm [41]; for instance, it includes implementations of 
(variations of) Z's "slow algorithm", and algorithms by F. Chyzak (and B. Salvy), and N. Takayama. 
In this context we will have to exploit closure properties of classes of special (resp. holonomic) 
sequences and functions; an introduction to computer algebra methods for the univariate case can 
be found in [15]. 

Moreover, the zb.m package strongly suggests that there are, in fact, four linearly independent 
virial recurrence relations, see more details on the article's website, but only three of them (e. g., 
(2.12)-(2.14)) are available in the literature. Another next challenge is to study the off-diagonal 
matrix elements that are important in applications [23], [24], [25], [29], [30], and [33]. 

Acknowledgment. We thank Doron Zeilberger for valuable discussions and encouragement. We 
are grateful to Christian Krattenthaler for his KTj^X style file for presenting Mathematica notebooks 
in this paper. 
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